Multi-objective optimization method for power system peak-shaving based on index linkage analysis

ABSTRACT

The invention relates to a multi-objective optimization method for power system peak-shaving based on index linkage analysis. The invention comprehensively considers five peak shaving indexes, and introduces a linkage analysis theory in economics to test stationarity and causality of data series of peak-shaving indexes. The purpose is to effectively clarify linkage relationship between indexes. A regression analysis method is used to quantify the mutual influence relationship. A multi-index coupling peak shaving optimization model for generation scheduling is developed and finally solved by using genetic algorithm. The invention can objectively describe comprehensive peak shaving requirements of a provincial power grid, and is helpful for generation dispatcher to reduce subjective assumptions in the decision-making. Compared with a conventional single-index peak-shaving model, peak-shaving effect of our method is obviously improved, which provides a novel technical way for peak-shaving generation scheduling and decision-making in short-term operations of power systems.

TECHNICAL FIELD

The invention relates to the field of hydropower scheduling, focusing on a multi-objective optimization method for power system peak-shaving based on index linkage analysis.

BACKGROUND

Peak-shaving is a core task of the short-time generation scheduling of power systems, and it's also one of the key issues at home and abroad. Due to the large load difference between peak and off-peak and the increasing intermittent energy grid connection, most provincial power grids in China are facing severe peak-shaving pressure. Peak-shaving has become one of the main targets in the daily power generation scheduling of these power grids. From the perspective of the power grid, the main purpose of peak-shaving is optimizing the power generation plans of the dispatchable power plants to reduce the peak-valley difference of residual loads and smooth the fluctuations in the net load curve for coal-fired units with poor regulating performance. In order to achieve this goal, how to accurately describe the peak-shaving requirements becomes a key part of power generation scheduling.

Nowadays, the technical measures to deal with the peak-shaving pressure of power systems mainly include two aspects: One is to focus on the construction of energy storage systems and its potential of peak-shaving, such as introducing a peak demand charge into the energy cost function, the joint operation optimization of a hybrid system including a photovoltaic farm and battery-based energy storage, the complementary operation of photovoltaic-wind plant and hydro-pumping storage. The other is the peak-shaving optimization of conventional power plants such as thermal power plants, hydropower plants and so on. This kind of methods generally formulates models based on the variable of the net load. Some forms of a function of the net load are used to describe the peak-shaving needs, such as the minimum of the maximum residual peak-valley difference, the minimum of the variance or mean square deviation of remaining load, the maximum of the peak-shaving capacity and so on. Overall, these modelling methods can effectively reduce the maximum peak-valley difference or smooth the fluctuations in the load curve in appropriate situations. However, for the peak-shaving pressure or requirements facing provincial power grids, it is usually described by a peak-shaving index such as the peak-valley difference, the mean square deviation of remaining load, etc. The comprehensive analysis of the load demand is lacked. When the load demand of different power girds changed at different dates, there is a lack of attention on the applicability of each peak-shaving index and which peak-shaving index should be adopted to optimize the generation schedules for better peak-shaving results. It's a very important problem in the research of peak operations.

For the above problem, the invention proposes a multi-objective optimization method for power system peak-shaving based on index linkage analysis. This work is supported by the National Natural Science Foundation of China (52079014). Based on the engineering background of Xiluodu hydropower plant transmitting power to ZheJiang Power Grid, the invention is tested and applied. The results show that the invention can effectively avoid the deviation of peak-shaving results caused by single-index optimization as well as the uncertainty and irrationality of subjective factors. Thus, the optimized generation schedules are more effective and practical.

SUMMARY

The technical problem to be solved by the invention is to provide a multi-objective optimization method for power system peak-shaving based on index linkage analysis. The purpose is to reduce uncertainty and irrationality caused by the weight coefficient selection in the multi-objective processing strategy and to obtain balanced peak-shaving optimization results.

The technical solution of the invention:

A multi-objective optimization method for power system peak-shaving based on index linkage analysis is characterized by the following steps:

(1) Set initial calculation conditions, including operation conditions and constraints of hydropower plants, as well as conditions for electric and hydraulic operation demands;

(2) Five peak-shaving indexes are introduced, including load variance

${F_{1} = {\sum\limits_{i = 1}^{n}\left( {P_{r,i} - {\overset{¯}{P}}_{r}} \right)^{2}}},$

load difference between peak and off-peak F₂=P_(r,max)−P_(r,min), maximum variation rate of loads F₃=|P_(r,i)−P_(r,i−1)|_(max), maximum load F₄=P_(r,max), minimum load F₅=P_(r,min); where P_(r) is residual load series; P_(r,i) is residual load in period i; P_(r,max) is maximum residual load; P_(r,min) is minimum residual load; n is total number of periods. Thus, a multi-objective function is obtained:

$\begin{matrix} {{MOP}\left\{ \begin{matrix} {\min F}_{1} \\ {\min F}_{2} \\ {\min F}_{3} \\ {\min F}_{4} \\ {\max F}_{5} \end{matrix} \right.} & (1) \end{matrix}$

(3) According to historical data of daily load curve, values of the above five indexes are calculated to obtain five time series that are named l₁, l₂, l₃, l₄ and l₅, respectively.

(4) The following formula is used to calculate correlation coefficients between pairs of five series mentioned above:

$\begin{matrix} {r = \frac{\sum\limits_{i = 1}^{n}{\left( {l_{x,i} - \overset{¯}{l_{x}}} \right)\left( {l_{y,i} - \overset{¯}{l_{y}}} \right)}}{\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{x,,i} - \overset{¯}{l_{x}}} \right)^{2}}\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{y,i} - \overset{¯}{l_{y}}} \right)^{2}}}} & (2) \end{matrix}$

-   where r is correlation coefficient; l_(x,i) and l_(y,i) are values     of sequence l_(x) and sequence l_(y) in period i, respectively;     l_(x) and l_(y) are average values of sequence l_(x) and sequence     l_(y), respectively. Correlation between five indexes is judged     according to the follows: two indexes are highly correlated if     0.8<|r|<1; two indexes are moderately correlated if 0.3<|r|<0.8 and     two indexes are lowly correlated if |r|<0.3.

(5) Based on the correlation between indexes, correlation analysis is carried out for indexes with moderate or higher correlation; a weight coefficient method is adopted for indexes with low correlation; a multi-objective peak-shaving objective is transformed into a single-objective peak-shaving model; specific steps are as follows:

Step1. Check stationarity of time series. ADF test method is used; a generation process of sequence l is one of the following three forms:

$\begin{matrix} {{{\Delta\; l_{t}} = {{\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (3) \\ {{{\Delta\; l_{t}} = {\mu + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (4) \\ {{{\Delta l_{t}} = {\mu + {\alpha t} + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (5) \end{matrix}$

where αt is time trend; μ is displacement item;

$\sum\limits_{i = 1}^{k}{\gamma_{i}{\Delta l}_{t - i}}$

is lag item; k is number of lag items; Δl_(t) is increment of t in sequence l that is calculated by Δl_(t)=l_(t)−l_(t−1); original hypothesis and alternative hypothesis are

$\left\{ {\begin{matrix} {{H_{0}:\beta} = {0\left( y_{t} \right.}} & \left. {nonstationary} \right) \\ {H_{1}:{\beta < {0\left( y_{t} \right.}}} & \left. {stationary} \right) \end{matrix},} \right.$

respectively; a judging criteria is

$\left\{ {\begin{matrix} {{{{ADF}\ {value}} \geq {threshold}},{{accept}\mspace{14mu} H_{0}},y_{t}} & {nonstationary} \\ {{{{ADF}\ {value}} < {threshold}},{{accept}\mspace{14mu} H_{0}},y_{t}} & {stationary} \end{matrix};} \right.$

test the above three forms in turn; if one of them rejects original hypothesis, the sequence would be stable;

Step 2. Granger causality test is conducted among sequences with moderate correlation and stability. Regression models between them are as follows:

$\begin{matrix} {l_{y,t} = {{\sum\limits_{i = 1}^{p}{\alpha_{l}l_{X/l}}} + {\sum\limits_{i = 1}^{p}{\beta_{i}l_{y,{t - i}}}} + ɛ_{1t}}} & (6) \\ {l_{x,t} = {{\sum\limits_{i = 1}^{p}{\lambda_{i}l_{y,{t - i}}}} + {\sum\limits_{i = 1}^{p}{\delta_{i}l_{x,{t - i}}}} + ɛ_{2t}}} & (7) \end{matrix}$

where α_(i), β_(i), γ_(i) and δ_(i) are regression coefficients; ε_(1t), ε_(2t) are random error items; p is a displacement item; carry out hypothesis testing for the two models, respectively; original hypothesis and alternative hypothesis are

$\begin{matrix} \left\{ \begin{matrix} {{H_{0}:\alpha_{i}} = 0} & \left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right) \\ {H_{1}:{\alpha_{i} \neq 0}} & \left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right) \end{matrix} \right. \\ {and} \\ \left\{ \begin{matrix} {{H_{0}:\lambda_{1}} = 0} & \left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right) \\ {H_{1}:{\lambda_{1} \neq 0}} & \left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right) \end{matrix} \right. \end{matrix}$

respectively; according to test results, the Granger causality between sequences could be obtained; if there is Granger causality between l_(x) and l_(y), then corresponding index X would has an impact on index Y.

Step 3. Co-integration of unstable sequences. If an unstable sequence l_(y) becomes a stable sequence after D times difference, then l_(y) is integrated of order D; use OLS co-integration regression equation for two unstable sequences l_(x) and l_(y) with same integrated order; where T, U are regression coefficients; X_(t) is residual; when X_(t) is stable, then corresponding index X would has an impact on index Y;

Step 4. Regression analysis based on Granger causality of stable sequence and quantitative variation relationship of unstable sequences; choose the most affected one of five indexes as θ₀; indexes which affect θ₀ are denoted as θ₁, θ₂, . . . , θ_(m), respectively; where m is number of indexes that affect θ₀; a formula θ₀=f (θ₁, θ₂, . . . , θ_(m)) is calculated by regression analysis;

Step 5. Construct a peak-shaving objective function; based on formula θ₀=f (θ₁, θ₂, . . . , θ_(m)), the multi-objective function mentioned above can be transformed into a single-objective function by using a coefficient weight method for residual indexes such as θ_(n+1), θ_(n+2), . . . , θ₄; the function expression is as follows:

$\begin{matrix} {{\min\; F} = {{f\left( {\theta_{1},\theta_{2},\ldots\mspace{14mu},\theta_{m}} \right)} + {\sum\limits_{i = {m + 1}}^{4}{\beta_{i}\theta_{i}}}}} & (8) \end{matrix}$

(5) Get the generation profile of power plants by using Genetic Algorithm method to solve the objective function in (4).

The invention has the following beneficial effects: compared with a normal single-objective peak-shaving model, a multi-objective coupling peak-shaving models could quantitatively describe reasonable relationship among five indexes of fluctuation degree, fluctuation range, peak value, valley value and maximum variation rate of load, balance peak-shaving results; it could effectively avoid an error of peak-shaving caused by a single-objective optimization; different indexes essentially reflect different peak-shaving demands of a power grid; through a correlation analysis between indexes, it could specify and definite load regulation demands of different power grids and different times, and obtain a more reasonable peak-shaving model, making optimized generation schedules more effective and useful; the invention describes complicated peak-shaving demand of power systems by using economic theory methods; this research approach is benefit to promote description of practical problems in power system generation scheduling, which could provide relatively objective optimization results and effectively avoid uncertainty and irrationality caused by subjective factors.

DESCRIPTION OF DRAWINGS

FIG. 1: A general solution frame diagram of the invention.

FIG. 2: A causal diagram of indexes.

FIG. 3: Analysis of fitting effect.

DETAILED DESCRIPTION

The specific embodiments of the present invention are further described below in conjunction with the drawings and technical solutions.

In order to describe peak-shaving effect accurately, the invention comprehensively considers variance, peak-valley difference, maximum, minimum and maximum variation rate of residual loads. Among them, the variance mainly reflects frequency of load changes in rise and fall directions, which is a requirement for flexible adjustment ability of generating units:

$\begin{matrix} {\sigma^{2} = {\sum\limits_{\iota - 1}^{n}\left( {P_{r,\iota} - {\overset{¯}{P}}_{r}} \right)^{2}}} & (9) \\ {P_{r,i} = {P_{i} - N_{i}}} & (10) \end{matrix}$

where P_(r,i) is residual load at period i; P_(i) is load of receiving-end power grid; N_(i) is generation of optimized power plant.

The peak-valley difference Δ is a variation range of load, the smaller the range is, the smaller the peak-shaving capacity demand of generating units is:

Δ=P _(r,max) −P _(r,min)  (11)

-   where P_(r,max) is maximum residual load; P_(r,min) is minimum     residual load.

The maximum load N_(r,max) reflects peak-shaving pressure of peak load periods that is a period with the maximum load response pressure of generating units. To ensure safe and stable operations of power grids, the residual load is expected to stay at a low level during peak periods in order to avoid peak-shaving difficulty of generating units.

$\begin{matrix} {N_{{ma}\; x} = P_{r,{{ma}\; x}}} & (12) \end{matrix}$

The minimum load N_(r,min) reflects peak-shaving pressure of valley load periods, which is limited by a forced generation of non-regulating units. During low load periods, the output of generating units is difficult to be effectively reduced with load, and it forces thermal units to work with in-depth peak-shaving, thereby resulting in a significant increase in energy consumption. Therefore, the residual load is expected to stay at a high level during valley periods.

$\begin{matrix} {N_{m\; i\; n} = P_{r,{m\; i\; n}}} & (13) \end{matrix}$

The maximum variation rate ΔP_(r,max) of load is maximum variation of load in adjacent periods that reflects demand of load for climbing ability of generating units. The smaller the maximum variation rate of load is, the lower the demand for climbing ability of generating units is.

$\begin{matrix} {{{\Delta\; P_{r,{m\;{ax}}}\mspace{11mu}\text{=}\mspace{11mu} P_{r,i}} - P_{r,{i - 1}}}}_{m\;{ax}} & (14) \end{matrix}$

The objective function is constructed to minimize the variance, peak-valley difference, maximum and maximum variation rate of residual loads, and maximize the minimum of residual loads to ensure peak-shaving demands of different loads.

$\begin{matrix} {{MOP}\left\{ \begin{matrix} {\min\;\sigma^{2}} \\ {\min\;\Delta} \\ {\min\; N_{m\;{ax}}} \\ {\max\; N_{m\; i\; n}} \\ {\min\;\Delta\; P_{r,{m\;{ax}}}} \end{matrix} \right.} & (15) \end{matrix}$

The objective function mentioned above should meet the following constraints:

Water balance equation:

$\begin{matrix} {V_{t + 1} = {V_{t} + {3600\left( {Q_{t} - S_{t}} \right)\Delta t}}} & (16) \end{matrix}$

where S_(t)=q_(t)+d_(t);

V_(t) is reservoir capacity of power plant at time t; Q_(t) is reservoir inflow in period t, m³/s; S_(t) is storage outflow of power plant in period t, m³/s; q_(t) is power discharge of power plant in period t, m³/s; d_(t) is abandoning water of power plant in period t, m³/s.

Relationship of water level and reservoir capacity:

$\begin{matrix} {Z_{t} = {f_{ZV}\left( V_{t} \right)}} & (17) \end{matrix}$

where Z_(t) is water level of power plant in period t, m; f_(ZV) (•) is a relation curve of water level and reservoir capacity of power plant.

The relationship of tailrace level and discharge capacity:

$\begin{matrix} {Z_{t}^{d} = {f_{z^{d}}\left( S_{t} \right)}} & (18) \end{matrix}$

where Z_(t) ^(d) is tailrace level of power plant in period t, m; f_(Z) _(d) (•) is a relation curve of tailrace level and discharge capacity of power plant.

Water head relationship:

$\begin{matrix} {h_{t} = {{\left( {Z_{t} + Z_{t + 1}} \right)/2} - Z_{t}^{d}}} & (19) \end{matrix}$

where h_(t) is water head of power plant in period t, m.

The relationship of power station output:

N _(t) =f _(qh) (q _(t) , h _(t))  (20)

where N_(t) is generation of power plant in period t, MW; f_(qh) (•) is a NQH curve of power plant.

Electricity generation control:

$\begin{matrix} {E_{total} = {\sum\limits_{t = 1}^{T}{N_{t} \times \Delta t}}} & (21) \end{matrix}$

where E_(total) is given total amount of power generation, MWh.

Constraint of reservoir water level:

$\begin{matrix} {{\underset{\_}{Z}}_{t} \leq Z_{t} \leq {\overset{\_}{Z}}_{t}} & (22) \end{matrix}$

where Z _(t) and Z _(t) are upper and lower water level of reservoir in period t, respectively, m.

Upper and lower limits of power generation:

$\begin{matrix} {{\underset{\_}{N}}_{t} \leq N_{t} \leq {\overset{\_}{N}}_{t}} & (23) \end{matrix}$

where N _(t) and N _(t) are upper and lower average generation in period t of power plant, respectively, MW.

Constraint of power discharge:

$\begin{matrix} {{\underset{¯}{q}}_{t} \leq q_{t} \leq {\overset{\_}{q}}_{t}} & (24) \end{matrix}$

where q _(t) and q _(t) are upper and lower power discharge in period t of power plant, respectively, m³/s.

Constraint of storage outflow:

$\begin{matrix} {{\underset{\_}{S}}_{t} \leq S_{t} \leq {\overset{\_}{S}}_{t}} & (25) \end{matrix}$

where S _(t) and S _(t) are upper and lower storage outflow in period t of m power plant, m³/s.

Taking the above multi-objective optimization model as the background, the invention is applied in a complete way according to the following (1) to (5) steps:

(1) Set initial calculation conditions, including operation conditions and constraints of hydropower plants, as well as conditions for electric and hydraulic operation demands;

(2) Five peak-shaving indexes are introduced, including load variance

${F_{1} = {\sum\limits_{i = 1}^{n}\left( {P_{r,i} - {\overset{¯}{P}}_{r}} \right)^{2}}},$

load difference between peak and off-peak F₂=P_(r,max)−P_(r,min), maximum variation rate of loads F₃=|P_(r,i)−P_(r,i−1)|_(max), maximum load F₄=P_(r,max), minimum load F₅=P_(r,min); where P_(r) is residual load series; P_(r,i) is residual load in period i; P_(r,max) is maximum residual load; P_(r,min) is minimum residual load; n is total number of periods. Thus, we get a multi-objective function:

$\begin{matrix} {{MOP}\left\{ \begin{matrix} {\min\; F_{1}} \\ {\min\; F_{2}} \\ {\min\; F_{3}} \\ {\min\; F_{4}} \\ {\max\; F_{5}} \end{matrix} \right.} & (26) \end{matrix}$

(3) According to historical data of daily load curve, we can calculate values of the above five indexes to obtain five time series that are named l₁, l₂, l₃, l₄ and l₅, respectively.

(4) The following formula is used to calculate correlation coefficients between pairs of five series mentioned above:

$\begin{matrix} {r = \frac{\sum\limits_{i = 1}^{n}{\left( {l_{x,i} - \overset{\_}{l_{x}}} \right)\left( {l_{y,i} - \overset{\_}{l_{y}}} \right)}}{\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{x,i} - \overset{\_}{l_{x}}} \right)^{2}}\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{y,i} - \overset{\_}{l_{y}}} \right)^{2}}}} & (27) \end{matrix}$

where r is correlation coefficient; l_(x,i) and l_(y,i) are values of sequence l_(x) and sequence l_(y) in period i, respectively; l_(x) and l_(y) are average values of sequence l_(x) and sequence l_(y), respectively. Correlation between five indexes is judged according to the follows: two indexes are highly correlated if 0.8<|r|<1; two indexes are moderately correlated if 0.3<|r|<0.8 and two indexes are lowly correlated if |r|<0.3.

(5) Based on the correlation between indexes, correlation analysis is carried out for indexes with moderate or higher correlation; a weight coefficient method is adopted for indexes with low correlation; a multi-objective peak-shaving objective is transformed into a single-objective peak-shaving model; specific steps are as follows:

Step1. Check stationarity of time series. ADF test method is used; a generation process of sequence l is one of the following three forms:

$\begin{matrix} {{{\Delta\; l_{t}} = {{\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (28) \\ {{{\Delta l_{t}} = {\mu + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (29) \\ {{{\Delta\; l_{t}} = {\mu + {\alpha t} + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (30) \end{matrix}$

where αt is time trend; μ is displacement item;

$\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}$

is lag item; k is number of lag items; Δl_(t) is increment of t in sequence l that is calculated by Δl_(t)=l_(t)−l_(t−1); original hypothesis and alternative hypothesis are

$\left\{ {\begin{matrix} {{H_{0}\text{:}\mspace{14mu}\beta} = {0\left( {y_{t}\mspace{14mu}{nonstationary}} \right)}} \\ {{H_{1}\text{:}\mspace{11mu}\beta} < {0\left( {y_{t}\mspace{14mu}{stationary}} \right)}} \end{matrix},} \right.$

respectively; a judging criteria is

$\left\{ {\begin{matrix} {{{{ADF}\mspace{14mu}{value}} \geq {threshold}},{{accept}\mspace{14mu} H_{0}},{y_{t}\mspace{14mu}{nonstationary}}} \\ {{{{ADF}\mspace{14mu}{value}} < {threshold}},{{accept}\mspace{14mu} H_{0}},{y_{t}\mspace{14mu}{stationary}}} \end{matrix};} \right.$

test the above three forms in turn; if one of them rejects original hypothesis, the sequence would be stable;

Step 2. Granger causality test is conducted among sequences with moderate correlation and stability. Regression models between them are as follows:

$\begin{matrix} {l_{y,t} = {{\sum\limits_{i = 1}^{p}{\alpha_{i}l_{x,{t - i}}}} + {\sum\limits_{i = 1}^{p}{\beta_{i}l_{y,{t - i}}}} + ɛ_{1t}}} & (31) \\ {l_{xt} = {{\sum\limits_{i = 1}^{p}{\lambda_{i}l_{y,{t - i}}}} + {\sum\limits_{i = 1}^{p}{\delta_{i}l_{x,{t - i}}}} + ɛ_{2t}}} & (32) \end{matrix}$

where α_(i), β_(i), γ_(i) and δ_(i) are regression coefficients; ε_(1t), ε_(2t) are random error items; p is a displacement item; carry out hypothesis testing for the two models, respectively; original hypothesis and alternative hypothesis are

$\left\{ {\begin{matrix} {{H_{0}:\mspace{14mu}\alpha}_{i} = {0\left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}} \\ {{H_{1}:\mspace{14mu}\alpha}_{i} \neq {0\left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}} \end{matrix}{and}\mspace{14mu}\left\{ \begin{matrix} {{H_{0}:\mspace{14mu}\lambda}_{i} = {0\left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{y}\mspace{14mu}{by}\mspace{14mu} l_{x}} \right)}} \\ {{H_{1}:\mspace{14mu}\lambda}_{i} \neq {0\left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{y}\mspace{14mu}{by}\mspace{14mu} l_{x}} \right)}} \end{matrix} \right.} \right.$

respectively; according to test results, the Granger causality between sequences could be obtained; if there is Granger causality between l_(x) and l_(y), then corresponding index X would has an impact on index Y.

Step 3. Co-integration of unstable sequences. If an unstable sequence l_(y) becomes a stable sequence after D times difference, then l_(y) is integrated of order D; use OLS co-integration regression equation for two unstable sequences l_(x) and l_(y) with same integrated order; where T, U are regression coefficients; X_(t) is residual; when X_(t) is stable, then corresponding index X would has an impact on index Y;

Step 4. Regression analysis based on Granger causality of stable sequence and quantitative variation relationship of unstable sequences; choose the most affected one of five indexes as θ₀; indexes which affect θ₀ are denoted as θ₁, θ₂, . . . , θ_(m), respectively; where m is number of indexes that affect θ₀; a formula θ₀=f (θ₁, θ₂, . . . , θ_(m)) is calculated by regression analysis;

Step 5. Construct a peak-shaving objective function; based on formula θ₀=f (θ₁, θ₂, . . . , θ_(m)), the multi-objective function mentioned above can be transformed into a single-objective function by using a coefficient weight method for residual indexes such as θ_(n+1), θ_(n+2), . . . , θ₄; the function expression is as follows:

$\begin{matrix} {{\min\mspace{11mu} F} = {{f\left( {\theta_{1},\theta_{2},\ldots\mspace{11mu},\theta_{m}} \right)} + {\sum\limits_{i = {m + 1}}^{4}{\beta_{i}\theta_{i}}}}} & (33) \end{matrix}$

(5) Get the generation profile of power plants by using Genetic Algorithm method to solve the objective function in (4).

Now, the presented model was verified by applying it to an actual engineering of Xiluodu Hydropower Plant. Xiluodu is one of main power suppliers of “Electricity Transmission Project from West to East”, and it's also the biggest hydropower plant on the Jinsha River. There are 9 identical hydropower units of 770 MW each installed on the left and right sides of the river. Among them, generating units of the plant on the left banks provide electricity to Zhejiang via the Binjin UHVDC line. The proportion of thermal power in the Zhejiang Power Grid is up to 95%, which lacks flexible generation capacity and exhibits large daily load differences between peak and off-peak periods, facing a high peak-shaving pressure. Therefore, it's necessary to make use of flexible regulation characteristics of hydropower to meet peak-shaving requirements. Typical daily load curves of summer and winter in Zhejiang province are selected for calculation. The calculation results of correlation coefficient between indexes are given in Table 2, among which the maximum variation rate is of low correlation with other indexes, thus only four indexes, variance, peak-valley difference, maximum, minimum, need to be analyzed jointly. Stability test results are listed in Table 3, and it can be obtained that four indexes are stable without co-integration. The results of Granger causality test are listed in Table 4, so a causal guidance diagram can be drawn between these indexes, as shown in FIG. 2. Regression analysis was conducted according to the causality, and the following fitting formula could be obtained, and the fitting effect is shown as FIGS. 3:

θ₀=9.23×10⁸−0.041×θ₁ ²+0.462×θ₂ ²−37535×θ₂+0.132×θ₃ ²−9529×θ₃

-   where θ₀ is variance; θ₁ is peak-valley difference; θ₂ is minimum;     θ₃ is maximum.

The index of maximum variation rate is treated by a weight coefficient method, and the following function could be obtained:

min F=9.23×10⁸−0.041×θ₁ ²+0.462×θ₂ ²−37535×θ₂+0.132×θ₃ ²−9529×θ₃+βθ₄

-   where θ₄ is maximum variation rate; β is weight coefficient.

Using the objective function mentioned above for the optimization, the indexes and decreasing ranges of corresponding index could be obtained and listed in Table 5 and Table 6, respectively. Among them, indexes of variance, peak-valley difference, maximum and the maximum variation rate showed a significant decline and the index of minimum showed a small decline. It can be seen that it has a significant peak-shaving effect in summer and winter, and effectively alleviates the peak-shaving pressure of the power grid.

TABLE 1 Basic data of power plant Normal Maximum Maximum Plant Dead water high water Installed Regulation reservoir power name level/m level/m capacity/MW performance discharge/m³/s discharge/m³/s Xiluodu 374.00 600.00 13860 multi-year 100000 7749 regulating

TABLE 2 Calculation results of correlation coefficient Maximum Fluctuation Fluctuation Peak Valley variation R degree range value value rate Variance 1 Peak-valley 0.949 1 difference Maximum 0.720 0.861 1 Minimum −0.950 −0.922 −0.597 1 Maximum 0.140 0.276 0.094 −0.056 1 variation rate

TABLE 3 Stability testing ADF test 1% critical 5% critical value value value C1onclusion Variance −4.42021 −3.43927 −2.8653 stable Peak-valley difference −3.74418 −3.43926 −2.86536 stable Maximum −4.4423 −3.43926 −2.86536 stable Minimum −4.13918 −3.43926 −2.86536 stable

TABLE 4 Results of Granger F- Null hypothesis statistic P value Peak-valley difference is not Granger cause 25.349 5  3.00 × 10−27 of variance Variance is not Granger cause of Peak- 13.853 4  6.00 × 10−15 valley difference Minimum is not Granger cause of maximum 2.335 3 0.030 6 Variance is not Granger cause of minimum 1.616 2 0.139 8 Maximum is not Granger cause of variance 18.092 4  1.00 × 10−19 Variance is not Granger cause of maximum 2.499 2 0.021 2 Maximum is not Granger cause of Peak- 4.781 8 9.00 × 10−05 valley difference Peak-valley difference is not Granger cause 0.488 7 0.817 1 of maximum Minimum is not Granger cause of Peak- 4.781 8 9.00 × 10−05 valley difference Peak-valley difference is not Granger cause 14.320 0  2.00 × 10−15 of minimum Minimum is not Granger cause of maximum 0.488 7 0.817 1 Maximum is not Granger cause of minimum 14.320 0  2.00 × 10−15

TABLE 5 Calculation results of indexes Index θ₀/MW² θ₁/MW θ₂/MW θ₃/MW θ₄/MW Typical Original 17 612 306  13 262  37 982 51 244 4 696 summer load days Residual 4 520 608 6 333 37 982 44 315 1 616 load Typical Original 2 531 061 5 156 43 401 48 557 4 183 winter load days Residual   144 675 1 103 41 789 42 892 1 103 load

TABLE 6 Decreasing range of indexes (%) Index θ₀ θ₁ θ₂ θ₃ θ₄ Typical summer days 74.3 52.2 0.0 13.5 65.6 Typical winter days 94.3 78.6 3.7 11.7 73.6 

1. A multi-objective optimization method for power system peak-shaving based on index linkage analysis, wherein comprising the following steps: (1) set initial calculation conditions, including operation conditions and constraints of hydropower plants, as well as conditions for electric and hydraulic operation demands; (2) five peak-shaving indexes are introduced, including load variance ${F_{1} = {\sum\limits_{i = 1}^{n}\left( {P_{r,i} - {\overset{¯}{P}}_{r}} \right)^{2}}},$ load difference between peak and off-peak F₂=P_(r,max)−P_(r,min), maximum variation rate of loads F₃=|P_(r,i)−P_(r,i−1)|_(max), maximum load F₄=P_(r,max), minimum load F₅=P_(r,min); where P_(r) is residual load series; P_(r,i) is residual load in period i; P_(r,max) is maximum residual load; P_(r,min) is minimum residual load; n is total number of periods; thus, a multi-objective function is obtained: $\begin{matrix} {{MOP}\left\{ \begin{matrix} {\min F}_{1} \\ {\min F}_{2} \\ {\min F}_{3} \\ {\min F}_{4} \\ {\max F}_{5} \end{matrix} \right.} & (1) \end{matrix}$ (3) according to historical data of daily load curve, values of the above five indexes are calculated to obtain five time series that are named l₁, l₂, l₃, l₄ and l₅, respectively; (4) Eq. (2) is used to calculate correlation coefficients between pairs of five series mentioned above: $\begin{matrix} {r = \frac{\sum\limits_{i = 1}^{n}{\left( {l_{x,i} - {\overset{\_}{l}}_{x}} \right)\left( {l_{y,i} - {\overset{\_}{l}}_{y}} \right)}}{\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{x,i} - {\overset{\_}{l}}_{x}} \right)^{2}}\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{y,i} - {\overset{\_}{l}}_{y}} \right)^{2}}}} & (2) \end{matrix}$ where r is correlation coefficient; l_(x,i) and l_(y,i) are values of sequence l_(x) and sequence l_(y) in period i, respectively; l_(x) and l_(y) are average values of sequence l_(x) and sequence l_(y), respectively; correlation between five indexes is judged according to the follows: two indexes are highly correlated if 0.8<|r|<1; two indexes are moderately correlated if 0.3<|r|<0.8 and two indexes are lowly correlated if |r|<0.3; (5) based on the correlation between indexes, correlation analysis is carried out for indexes with moderate or higher correlation; a weight coefficient method is adopted for indexes with low correlation; a multi-objective peak-shaving objective is transformed into a single-objective peak-shaving model; specific steps are as follows: Step1. check stationarity of time series; ADF test method is used; a generation process of sequence l is one of the following three forms: $\begin{matrix} {{{\Delta l_{t}} = {{\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (3) \\ {{{\Delta l_{t}} = {\mu + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (4) \\ {{{\Delta l_{t}} = {\mu + {\alpha t} + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (5) \end{matrix}$ where αt is time trend; μ is displacement item; $\sum\limits_{i = 1}^{k}{\gamma_{i}{\Delta l}_{t - i}}$ is lag item; k is number of lag items; Δl_(t) is increment of t in sequence l that is calculated by Δl_(t)=l_(t)−l_(t−1); original hypothesis and alternative hypothesis are $\left\{ {\begin{matrix} {{H_{0}:\mspace{14mu}\beta} = {0\left( {y_{t}\mspace{14mu}{nonstationary}} \right)}} \\ {H_{1}:\mspace{14mu}{\beta < \left( {y_{t}\mspace{14mu}{stationary}} \right)}} \end{matrix},} \right.$ respectively; a judging criteria is $\left\{ {\begin{matrix} {{{{ADF}\mspace{14mu}{value}} \geq {threshold}},{{accept}\mspace{14mu} H_{0}},{y_{t}\mspace{14mu}{nonstationary}}} \\ {{{{ADF}\mspace{14mu}{value}} < {threshold}},{{accept}\mspace{14mu} H_{0}},{y_{t}\mspace{14mu}{stationary}}} \end{matrix};} \right.$ test the above three forms in turn; if one of them rejects original hypothesis, the sequence would be stable; Step
 2. Granger causality test is conducted among sequences with moderate correlation and stability; regression models between them are as follows: $\begin{matrix} {l_{y,t} = {{\sum\limits_{i = 1}^{p}{\alpha_{i}l_{x,{t - i}}}} + {\sum\limits_{i = 1}^{p}{\beta_{i}l_{y,{t - i}}}} + ɛ_{1t}}} & (6) \\ {l_{x,t} = {{\sum\limits_{i = 1}^{p}{\lambda_{i}l_{y,{t - i}}}} + {\sum\limits_{i = 1}^{p}{\delta_{i}l_{x,{t - i}}}} + ɛ_{2t}}} & (7) \end{matrix}$ where α_(i), β_(i), γ_(i) and δ_(i) are regression coefficients; ε_(1t), ε_(2t) are random error items; p is a displacement item; carry out hypothesis testing for the two models, respectively; original hypothesis and alternative hypothesis are $\left\{ {\begin{matrix} {{H_{0}:\mspace{14mu}\alpha}_{i} = {0\left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}} \\ {{H_{1}:\mspace{14mu}\alpha}_{i} \neq {0\left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}} \end{matrix}{and}\mspace{14mu}\left\{ \begin{matrix} {{H_{0}:\mspace{14mu}\lambda}_{i} = {0\left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{y}\mspace{14mu}{by}\mspace{14mu} l_{x}} \right)}} \\ {{H_{1}:\mspace{14mu}\lambda}_{i} \neq {0\left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{y}\mspace{14mu}{by}\mspace{14mu} l_{x}} \right)}} \end{matrix} \right.} \right.$ respectively; according to test results, the Granger causality between sequences could be obtained; if there is Granger causality between l_(x) and l_(y), then corresponding index X would has an impact on index Y; Step
 3. co-integration of unstable sequences; if an unstable sequence l_(y) becomes a stable sequence after D times difference, then l_(y) is integrated of order D; use OLS co-integration regression equation for two unstable sequences l_(x) and l_(y) with same integrated order; where T, U are regression coefficients; X_(t) is residual; when X_(t) is stable, then corresponding index X would has an impact on index Y; Step
 4. regression analysis based on Granger causality of stable sequence and quantitative variation relationship of unstable sequences; choose the most affected one of five indexes as θ₀; indexes which affect θ₀ are denoted as θ₁, θ₂, . . . , θ_(m), respectively; where m is number of indexes that affect θ₀; a formula θ₀=f (θ₁, θ₂, . . . , θ_(m)) is calculated by regression analysis; Step
 5. construct a peak-shaving objective function; based on Eq. θ₀=f (θ₁, θ₂, . . . , θ_(m)), the multi-objective function mentioned above can be transformed into a single-objective function by using a coefficient weight method for residual indexes such as θ_(n+1), θ_(n+2), . . . , θ₄; it is shown as follows: $\begin{matrix} {{\min F} = {{f\left( {\theta_{1},\theta_{2},{.\;.\;.}\;,\;\theta_{m}} \right)} + {\sum\limits_{i = {m + 1}}^{4}{\beta_{i}\theta_{i}}}}} & (8) \end{matrix}$ (5) get the generation profile of power plants by using genetic algorithm method to solve the objective function in (4). 